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I. INTRODUCTION AND CONTEXT. 


Many aspects of the behaviors of physical, chemical and biological systems can be nnder- 
stood simply in terms of the dynamics of the averaged state variables and their deterministic 
evolntion eqnations. Since snch systems typically involve very large nnmbers of ensemble 
members, or long time averaging, fluctnations are highly snppressed with respect to the 
mean behavior. However, for some cases random flnctnations do not simply add negligible 
noise to the averaged dynamics, instead they give rise to fnndamentally different behav¬ 
iors. The dynamics of the averaged variables are thus insufficient to capture the system’s 
behavior in the stochasticity-dominated regime. Classic examples of systems in this regime 
include critical phenomena in physics, genetic drift and extinction in biology and diffusion 
dominated reactions in physical chemistry 

It is increasingly appreciated that for biological systems at the cellular and molecular 
scale, fluctuations can lead to single-cell and single molecule behaviors that considerably 
deviate from naive ensemble-averaged expectations. This is because the underlying bio¬ 
chemical and biophysical processes often involve reactions with small numbers of reactants. 
Thus the inherently probabilistic nature of these processes cannot be ignored [T3HT9]. In 
addition to this “intrinsic” stochasticity, cells may also have additional sources of cell-to-cell 
variability, known as the “extrinsic noise” [ 20 ] • Qualitatively distinct behaviors may emerge 
when such stochastic fluctuations dominate the system dynamics; familiar biological exam¬ 
ples on the cellular scale include stochastic switching between different phenotypes EZH and 
stochastic resonances in neurobiology I221- 

Often sharp changes in cellular behavior are triggered by thresholded events, i.e., by 
the attainment of a threshold value of a relevant cellular or molecular dynamical variable. 
Since the governing variable itself typically undergoes noisy or stochastic dynamics, there 
is a corresponding variability in the times when the same change occurs in each cell of a 
population. This time is called the “hrst passage” time and the corresponding process is 
a “hrst passage” (FP) process, referring to the event when a random variable first passes 
the threshold value. Even seemingly simple processes, such as the transport of molecules 
through channels or multivalent binding also fall under the umbrella of the First Passage 
processes. 

While stochastic effects in copy number huctuations have received considerable attention 
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in recent years, both experimentally and theoretically dMa [23H29], the stochasticity in 
the outcomes and the corresponding noise in the timing of cellular and molecular events has 
not received comparable attention. In part, this is due to the experimental challenges in 
obtaining high quality time series data amenable to analysis for timing noise, which requires 
making in vivo measurements at the single cell level [301 EH. However, increasingly, this 
challenge is being overcome through rapid development of single-cell technologies which fa¬ 
cilitate making such observations [301 - ET] . On the molecular scale advances in measurement 
techniques are starting to provide direct insights into the single molecule transport, interac¬ 
tions and signalling processes on the nanoscale [12H1H]. These technological developments 
have made it apposite to now develop the FP formalism specihcally for addressing current 
problems in cellular and molecular biology, i.e., for establishing quantitative relations be¬ 
tween the timing noise in stochastic events and the corresponding underlying stochastic 
dynamics of the thresholded variables. 


Mathematical techniques for modeling and analyzing First Passage processes were pio¬ 
neered a few decades ago, in the context of non-equilibrium physical chemistry and chemical 
physics 0 [HHII]- Detailed descriptions can be found in several textbooks and reviews 
[HEliaElSnHSI]. These techniques are now increasingly being adapted to problems in cel¬ 
lular, molecular and population biology. The renewed interest in FP problems in biological 
contexts has been reflected in several new works summarizing various aspects of the applica¬ 
tions of FP theory to these problems [71 1521 - IM] . However, many fundamental and practically 
useful results remain scattered across the literature in somewhat disparate communities. 


In this review we hrst present and elucidate fundamentals of the FP formalism within 
a unihed conceptual framework, which naturally integrates the existing techniques. We 
then discuss applications thereof, with emphasis on the practical use of FP techniques in 
biophysical systems. Our focus here is on covering a diverse set of analytical techniques; 
the number of reviewed biological applications is thus limited, out of necessity. We focus 
on three specihc areas: channel transport; receptor binding and adhesion; and single-cell 
growth and division. 
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II. FRAMEWORK. 


The presentation in this section is partially drawn from these textbooks and reviews: 

[21111E1171E21E2]. 

A. Stochastic processes. 

We first review fundamentals of the theory of stochastic processes. The system dynamics 
are specified by the set of its states, {S'}, and the transitions between them, S —)■ S', where 
S, S' G (S). For example, the state S can denote the position of a Brownian particle, the 
numbers of molecules of different chemical species, or any other variable that characterizes 
the state of the system of interest. Here we restrict ourselves to processes for which the 
transition rates depend only on the system’s instantaneous state, and not the entirety of its 
history. Such memoryless processes are known as Markovian and are applicable to a wide 
range of systems. We also assume that the transition rates do not explicitly depend on time, 
a condition known as stationarity. In this review we make the standard assumption that the 
transitions between the states are Poisson distributed random processes. In other words, 
the probability of transitioning from state S" to state S in an infinitesimal interval, dt, is 
a{S, S')dt, where a{S, S') is the transition rate. 

Examples. For a Brownian particle diffusing along a line, the state S is defined by the 
particle position; the transition rate is 2D/d‘^, where D is the diffusion coefficient and d is 
the step length. For a set of radioactive atoms undergoing decay with rate k per atom, the 
state S is defined by the number, n, of atoms that have not decayed yet, and the transition 
rate from state n to state n — 1 is nn. For a system with N reacting chemical species, the 
system state is defined by the concentrations of each reactant, (xi... xn), and the transition 
rates are functions of these concentrations. 

B. Time evolution equation(s) for the system. 

We now summarize the equations that govern the dynamical evolution of the probability 
that the system is in state S at time f, which we denote by P(S', t). Typically, such equations 
are written in one of three formalisms: the Master Equation, the Fokker-Planck Equation, 
or the Stochastic Differential Equation; each is summarized below in turn. Details of the 
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derivations can be found in [2H1]. 


1. The Master Equation. 


The Master Equation (ME) is the most general of the three formalisms and comprises of 
a set of linear ordinary differential equations. The ME is derived as follows. The probability 
of being in state S' at a time t + dt, P{S, t + dt), is the sum of the following two terms. First, 
the probability, P{S,t), that the system was already in the state S at time t and remained 
there during dt. Second, the probability that the system was originally in some other state 
S' at time t, times the probability that the system transitioned from S' to S during dt. 
Combining these two terms one obtains 


P{S,t + dt) = P{S,t) 


1 -^a(S',S)dt 

S' 


^P(S',t)a(S, S')dt. 

S' 


Taking the limit dt —)■ 0, we get the forward Master Equation (FME), or simply the Master 
Equation (ME): 


dtP(S, t) = J2 “(S. S')P(S', t)-Y, a(S'. S)P{S, t). (1) 

S' S' 

The hrst term in Eq. is the probability flux into the state S while the second term is the 
flux out of S. The Master Equation can be compactly written in operator notation as 


d,P{t)=MfP{t), (2) 

where P{t) is a vector with the components P{S,t) and is a linear operator with the 
components 

{Mf)ss> = a{S, S') - 5ss' «(S", S). (3) 

S" 

Note that Eq. ([^ conserves probability, dt iJ2s ~ which is guaranteed since 

ss' ~ S') — q;(S", S') = 0. For systems with discrete states, the 

operator Aif is simply a matrix with the above components. 

Often physical systems have states that are either characterized by a continuous variable 
s, or can be conveniently viewed as a continuous limit of the discrete states S. In this 
case, the discrete probabilities, P{S,t), are replaced by the probability density, p{s,t), that 
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specifies the probability that the state of the system lies in an inhnitesimal region [s, s + ds] 
in s-space: P{[s,s + ds],t) = p{s,t)ds. The sums in the Master Equation Q are then 
replaced by the corresponding integrals: 

dtp{s,t) = j a{s, s')p{s',t)ds' — p{s,t) j a{s',s)ds'. (4) 

Both in the discrete and the continuous cases the formal solution to the Master Equation 
can be written in terms of the initial probability distribution, P{to), at the initial time, to, 
as 


P{t) = (5) 

For discrete state variables, Eq. (|^ simply requires exponentiation of a matrix, whereas for 
the case of continuous state variables, it generally requires solution of the integral equation 

0 . 

Examples. In the previously mentioned example of radioactively decaying atoms, the 
Master Equation is dtP{n,t) = {n + l)KP{n + l,f) — nKP{n,t), where k is the decay rate 
per atom. For a particle performing an unbiased random walk with jump length a and total 
jump rate r, the Master Equation is dtP{x, t) = ^P{x + a,t) + ^P{x — a, t) — rP{x, t). 

2. Kramers-Moyal expansion and the Fokker-Planek equation. 

In some cases where the state space is continuous, there is a sense of locality, and one 
can dehne a “distance” between two states s and s', 6 = s' — s. A familiar example of 
this scenario is a Brownian particle on a line, its instantaneous state being specihed by its 
coordinate. The Master Equation, Eq. Q, then can be rewritten as 

dtp{s,t) = j r{5, s-\-5)p{s + 5,t)d5 — p{s,t) J r{S,s)d6, (6) 

where r{6, s) = q:(s + 6 , s) is the rate of jumping over a distance 6 , away from the state s. 
If r{6, s) rapidly decays with increasing S, over the lengthscale of typical variations of the 
probability density, p{s,t), one can expand r{S,s + 5) and p{s + 6,t) around <5 = 0. Thus 
p{s + 6, t) = p{s, t) +6 dsp{s, t) +6‘^dgp{s, t)/2 + 0{6^) and r{6, s + 5) = r{6, s) + 5 dsr{6, s) + 
d‘^dgr{6, s)/2 + 0{5^)-, this is known as the Kramers-Moyal expansion. Long-tailed transition 
rates, r{5,s), result in anomalous diffusion, not addressed in this review [55]. Substituting 
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these expansions in Eq. ([^ and keeping terms till the second order in S, the minimnm order 
necessary for obtaining non-trivial diffnsion-like motion, we arrive at the following partial 
differential eqnation, known in physics literatnre as the Fokker-Planck Equation: 

dtp(s, t) = -ds {A{s)p{s, t)) + {B{s)p{s, t)) , (7) 

where the functions A(s) and B{s) are, respectively, the hrst and the second moments of 
the transition rate r((5, s): 

/ OO f*00 

r{6,s)5d5, B{s) = / r{5,s)5'^d5. (8) 

OO j —OO 

As we have noted previously, the total probability is conserved by the Master Equation. 
Analogousy, the Fokker-Planck equation, Eq. Q, conserves probability and can be written 
as a local continuity equation for the probability density, 

dtp{s,t) + dsJ{s,t) =0, ( 9 ) 

where the quantity J{s,t) = A{s)p{s,t) — ds{B{s)p{s,t))/2 is the probability current It is 
important to emphasize that the Fokker-Planck equation is an uncontrolled approximation 
to the full Master Equation, and can lead to different results mi. 

Physical interpretation. When the Fokker-Planck equation is used to describe the move¬ 
ment of a physical particle under the action of a force f{s), the equilibrium probability 
density has to satisfy the Boltzmann-Gibbs distribution, p{s) oc exp{— f{x)dx/kBT). 
Rewriting the probability current in the Fokker-Planck equation as 

J{s,t) = A{s)p{s,t) - ^B{s)dsp{s,t), (10) 

with A(s) = A(s) — ^dsB{s), we see that the equilibrium solution to the EPF, which should 
satisfy the condition J{s,t) = 0, is p{s,t) oc exp(f'^2A(x)/B(x)dx). Comparing with the 
Boltzmann-Gibbs distribution, this imposes the constraint A(s) = —B{s)f{s)/2kBT, which 
is known as the Einstein relation. Written this way, the current J (s, t) has a simple physical 
interpretation: the hrst term in the current is the drift, arising due to the action of the 
force and characterized by a velocity A(s), and the second term represents the diffusive hux 
(Pick’s law) with the diffusion coefficient D{x) = B{x)/2. 
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3. Langevin and the Stochastic Differential Equations. 


When the state variable is continuous, the stochastic evolution of the system can be 
thought of as deterministic motion with added random fluctuations. A familiar example of 
this case is the Langevin equation that describes the motion of a diffusing Brownian particle, 

xit) = ( 11 ) 

where f{x) is the deterministic force acting on the particle, ^ is a random force that mimics 
the effects of random jumps, and /i is the mobility. Note that x{t) is now a random variable. 
The formal connection between this representation and the probability density of the previ¬ 
ous section is provided by the relation p{s,t) = {S{x{t) — s))^, where the average is over all 
the realizations of the random force, 

It can be shown that with the choice of /i = D/kT and — t'), this 

equation is mathematically equivalent to the following Fokker-Planck equation, 

dtp{s,t) = -dsipf{s)p{s,t)) + Dd^p{s,t). (12) 


More generally, any Fokker-Planck equation of the form 

dtp{x, t) = -d:, {A{x)p{x, t)) + {B{x)p{x, t)) 

has an equivalent stochastic differential equation (SDE) of the form 

x{t) = Ao(x) Bo{x)x{t), 


(13) 


(14) 


with delta-correlated random term (x(^)x(fO) = ~ ^0- However, due to mathematically 

pathological properties of the function x(t) (it is nowhere differentiable), when Bo{x) de¬ 


pends on X, the Eq. (14) is not unambiguosly dehned. In general, its interpretation requires 


re-dehnition of the rules of differentiation and integration, and many different SDEs can be 
chosen to correspond to the same FPE, depending on the interpretation. 

Historically, the two major interpretations are from Ito and Stratonovich. In both 
these formulations, Bo{x) = B{x)^. However, Aq{x) = A{x) in Ito interpretation while 
Ao{x) = A[x) — \dxB{x) in the Stratonovich interpretation. From the practical perspec¬ 
tive, Ito interpretation allows one to simulate the SDE using the usual forward Euler scheme. 
However, special differentiation and integration rules are required for analytical calculations. 




On the other hand, Stratonovich interpretation allows using the regular rules of calculus but 
has to be simulated using implicit schemes. We emphasize that the Fokker-Planck equa¬ 
tion does not suffer from such ambiguity of interpretation; SDEs corresponding to different 
interpretations of the same Fokker-Planck equation lead to the same physical results PE!- 

C. Backward evolution equations. 

1. The Backward Master Equation. 

The general form of the Master Equation derived in Eq. Q is also known as the forward 
Master Equation (FME), since it describes the evolution from an initial state to a state at a 
later time. The Master Equation is linear in the probabilities P{S,t) (see Eq. [^. Thus, its 
solution with any general initial condition, P{S, to), can be obtained as a linear combination 
of the conditional probabilities P(S,tjSi,to), which are solutions to the Master Equation 
for the special initial conditions, P{S,to) = ds,Si- Mathematically, the P{S,t\Si,to) are the 
Green’s functions of the Master equation. 

The time evolution of these conditional probabilities, P{S,t\Si,to), can also be described 
by an alternative linear equation instead of the Forward Master Equation, known as the 
Backward Master Equation (BME), which is especially useful in the context of First Passage 
problems. The key to deriving the BME equation is to consider the first step out of the 
initial state St at time to, rather than the last step of the trajectory, leading to the state S, at 
time t. Similar to the derivation of the FME, the conditional probability, P{S,t\Si,to), can 
be written down as the sum of the probabilities of two mutually exclusive events: (i) that the 
system transitioned to a different state S' during the time interval dt, with the probability 
a{S', Si)dt, and then evolved to a state S by time t, with the probability P{S, t\S', to + dt), or 
(ii) that the system was still in state Si at time to+dt, with the probability 1—Xls' C({,S', Si)dt, 
and then by time t evolved to the state S, with the probability P{S, t\S', to + dt). Together, 
these terms yield 




(15) 


S' 


9 


The stationarity condition, i.e. the the lack of explicit dependence of the transition rates on 
time, implies that the conditional probability P{S, t\Si, to) is a function only of t — to- Thus, 
P{S, t|S", to + dt) = P(S', t — dt\S\ to) in the above equation and so 

P{S, t|5„ to) = S ^ 

+ X P{S, t - dt\S', to). (16) 

S' 

Taking the limit dt —)■ 0, we get 

dtP{S,t\Si, to) = a{S', Si)P{S, t\S', to) - P{S, t\Si, to) a(5', S,) 

S' S' 

( 17 ) 

S' 

The resulting equation is known as the Backward Master Equation (BME); Aif, is the 
Backward Master Operator with the components (Alb)^^, = a{S',S) — a(S", S'), 

and it is a transpose of the the forward operator Aif dehned in Eq. ([^. The Backward 
Master Equation can be also written in the operator form; 

dtP^it) = Mb ■ (18) 

where is a vector whose t-th component is P{S,t\Si): 

= {...P{S,t\Si),P{S,t\S2)...P{S,t\S^)...). 


2. The Backward Fokker-Planck equation. 


The Backward Master Equation can be extended to the continuous case, similar to the 


procedure applied to the FME in Section |II B 2| , and can be approximated by the corre¬ 
sponding Backward Fokker-Planck equation. 


dtP{x,t\xi,to) = A{xi)d^P{x,t\xi,to) + 7;B{xi)dlP{x,t\xuto)- 

/ 


(19) 


It is analogous to the forward Fokker-Planck equation, Eq. ([^, except that the differential 
operators on the right hand side act on the initial state, x*, instead of the current state, 
X, resulting in B{xi) being outside of the derivative sign. Note that the Backward Fokker- 
Planck equation conserves the overall probability; however, it cannot be written as a local 
continuity equation with respect to the initial position Xj. 
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D. First Passage Processes. 


We have now set up the framework required to address the First Passage (FP) problem, 
which can be stated as the following question. For a stochastic Markov process that starts 
from the initial state Si at time to, what is the distribution of times, r = t — to, at which 
the system arrives at the specihc state Sf for the first timel We denote the probability 
density of this First Passage Time distribution by F(r; S'/|S'j); due to stationarity it does 
not depend explicitly on to- 

Naively, one may be tempted to guess that the hrst passage time distribution should be 
proportional to the probability to be in state S'/ at time t, P{Sf,t\Si,to). However, this is 
incorrect because P{Sf, t\Si, to) contains contributions from trajectories in which the system 
has already visited the hnal state S'/ at other instances between times t and to- In other 
words, P{Sf,t\Si,to) over-counts the number of hrst passage trajectories. See Fig for a 
graphic representation of the hrst passage time problem. 


t/} 

D 


Time 



FIG. 1. The First Passage Time problem. Distributions of times when a system, starting from 
an initial state, Xi, first visits specified threshold, x/, can be found by considering an auxiliary 


problem, with an absorbing boundary condition at x/. (See Section IID ) Shown here are 6 sample 
trajectories starting from Xi (orange); the first passage time for each trajectory is marked by the 
dotted vertical line at the intersection of the trajectory with the threshold, x/. In the auxiliary 
problem, the trajectories continuing beyond the first visitation event (shown in gray) are irrelevant 
and should not be counted. 
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1. First Passage Processes using Forward Master and Fokker-Planck equations. 

Counting trajectories that have not visited the final state previously is a combinatorially 
complex problem. It can be solved by considering an auxiliary version of the original prob¬ 
lem, in which once the system arrives at the state Sf, it remains there indefinitely. Thus it is 
not allowed subsequent transitions to another state. In other words, one places an absorbing 
boundary condition at the state Sf, with transition rates out of Sf, a{S, Sf), being set equal 
to zero for all S. 

In this auxiliary problem we define the survival probability S(t, Sf\Si,to) as the proba¬ 
bility that the system has not yet been absorbed at 5*/ by time t, after starting from Si at 
t = to'. 


S{t,Sf\S„to)= PiS,t\S„to). (20) 

S^Sf 

Note that stationarity assumption dictates that iS is a function of t — to only. Since, by 
definition, the probability of reaching Sf in a time interval [to+r, to+r + dr] is F(r, Sf\Si)dT, 
the probability of reaching Sf by time t is F{t; Sf\Si)dT. In other words, the probability 
that the First Passage Time is larger than r is S(to + r, Sf\Si,to), and therefore S(to + 
T,Sf\Si,to) is the cumulative distribution of F{t, Sf\Si). Intuitively, it is clear that the 
survival probability S{t) decreases in time with the rate equal to the probability current 
into the absorbing state Sf: 

dtS(t, Sf\Si) = —J{Sf,t\Si). This provides a prescription for obtaining the First Passage 
Time distribution, F{Sf,T\Si), by solving the Forward Master Equation, which yields the 
probabilities P{S,t\Si), and hence the probability flux into the absorbing state. 

Formal derivation. These arguments can be put in a mathematically rigorous form. The 
survival probability is related to the FPT distribution as 

S{t,Sf\S,,to) = l- r '\{T-,Sf\Si)dT. (21) 

Jo 

Thus, the survival probability is the cumulative probability distribution for F(r; Sf\Si) and 

F{T-,Sf\S,) = -dtS{t,Sf\S„to)\t=to+r. (22) 

Using this with the Forward Master Equation, and keeping in mind that {J^f)s,Sf = 0 
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because a{S, Sf) = 0 for all S (see Eq. [^, 


dtS(t,Sf\Si,to) = ^ dtP{S,t\Si,tQ) = E E {Mf)s,s'P{S',t\Si,to) 

Sf^Sf Sj^Sf S'j^Sf 

S'^Sj S^Sj 


Y, P(S\t\S,,h) Y^M!)ss' - (Mf)s,.s- 

S'^Sf \ s 

- Y S')P{S',t\S,, to) = -J{Sf, t\S,, to). 

S' 


(23) 


We have used the facts that f)s,S' = 0 due to the conservation of probability and that 


(■^f)sfS' ~ *^0 section IIB 1). The quantity J in the last line is the probability 

current from all accessible states into Sf. Comparing with Eq. 22, this proves our heuristic 
assertion that F(r;S'/|S'j) = J(S'/, to +to)- 

This result can also be obtained for a continuous variable using the Fokker-Planck equa¬ 
tion (Eq. 0 ). shown below for a simple one-dimensional case. Putting to = 0 and assuming 
Xi < Xf (and thus p{x > Xf,t) = 0), 


F{t;sf\si) = -dtS{t,Sf\si) = -dt 


.Sf 


p{s, t\si)ds 


t=T 


dtp{s,t\si)ds] = / dsJ{s,t\si)ds = J{sf,t\si), 


(24) 


since the current at inhnity vanishes, i.e., J{oo,t\si) = 0. 

To summarize, in order to calculate the probability density of the First Passage Times to 
state Sf, from state Si, one needs to solve the Forward Master or Fokker-Planck equation for 
the auxiliary process with the absorbing boundary condition at S'/, obtain the probability 
current J{t) into the absorbing state Sf, which then provides the FPT distribution through 
the relation F(r; S'/|S'j) = J{Sf, to + T\Si, to). 


2. First Passage Processes using Backward Master and Fokker-Planck equations 

The First Passage Time distribution can also be calculated using the backward formalism 
of Eq. The crucial insight is that the survival probability, S, also satishes the Backward 
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Master Equation. Setting to = 0, 


d,S{t, S,\S.) = - Y, StP{S,t\Si) = -YY,<-^'>hs-P(S,t\S') 

Sj^Sf Sj^Sf S' 

= -^(X0s.s. Y = -E(-^‘)as''S((,S/|S'). (25) 

S' S^Sf S' 

where we have used the linearity of Aib- Taking another time derivative, we find that the 
FPT probability density also obeys the Backward Master Equation: 


s,\s,) = -a!s(t; s,|Si) = (xOs.s. -?(<; S/|S'). (26) 

S' 

Thus the First Passage Time distribution can be obtained by solving the Backward Master 
Equation (and correspondingly, the Backward Fokker-Planck equation). 

Although solving the Backward Master Equation is not necessarily easier than solving 
the Forward Master Equation, it provides a relatively easy way of calculating the moments 
of the F(f; Sf\Si). For instance, the Mean First Passage Time (MFPT), defined as 

poo 

TiSf\Si)= dTTF{T;Sf\S,), (27) 

Jo 

can be calculated by applying the Backward Master operator to both sides: 


E(-^'>)ss-r(S/|s') 

S' 



dTTY(MOs.S‘F{r;Sj\S') 

S' 



dr TdtF{T;Sf\S,) 



dTdrS{T, Sf\Si, 0) 


= S{oo,Sf\S,)-S{0,Sf\Si) = -l, 


(28) 


where we have used the fact that 5(0, S'/1S',) = 1 and 5(oo, Sf\Si) =0. In matrix form. 


Mb-T = -1, 


(29) 


where T is the vector whose Tth component is T(5/|Sj). 


Higher moments can be obtained by sequential application of the reasoning of Eq. (28). 
This obviates the need for solving the full time-dependent differential Master Equation. 
Instead one can simply find the solution to the much simpler set of linear algebraic equations 
satisfied by T(S'/|Sj). 
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For continuous variables, the MFPT obeys the corresponding Backward Fokker-Planck 
equation [21 El EHl E3], 

= -1. (30) 


Heuristics. One can examine the simple logic behind the cumbersome mathematics of the 
Backward Equations with the following simple example. Consider a symmetric and homo¬ 
geneous random walker on a lattice, hopping with equal rates, r, to the left or right, thus 
changing its position, x, by ±a. We wish to hnd the Mean First Passage time, T{xo), of 
the random walker arriving at x = 0, starting from some position xq. Following the above 
“backward” arguments, any trajectory from xq to 0 can be decomposed into two mutually 
exclusive families of paths: one consisting of hrst jumping to the left, i.e., to xq — a and then 
going to X = 0 from there, and the other in which the random walker hrst jumps to the 
right, to Xq -I- a, and then proceeds to x = 0. For either family the hrst step, being a random 
Poisson process, takes a time 1/r on average. Thus the MFPT is found by averaging over 
these two families of equiprobable trajectories: 

T{xo) = ^ (^“ + ^(^0 “ ^ 


Rearranging terms. 


'^T{xo -a) + ^r(a:i, + a) - rT{xo) = -1, 


(32) 


in agreement with equation (28) above. Taking the limit a —?• 0, we recover the equation 


satished by T(xo) in terms of the continuous case Backward Fokker-Planck operator: 


D 


dxl 


T(xo) = -1. 


(33) 


3. First Passage Processes with multiple absorbing states. 

More complicated First Passage Processes can be addressed within the same framework. 
A question that arises frequently is the following. What is the FPT distribution for reaching 
the state Sf for the hrst time, without passing through a set of other states {Sf/} before 
that? The answer to this question can be obtained following a similar prescription as above, 
by considering an auxiliary problem with absorbing boundary conditions at S'/ as well as the 
states {S'//}. The main diherence from the previous case is that the probability of reaching 


15 




the final state Sf is not eqnal to one anymore: some trajectories get to one of the states 
{Sf/} first and should not be counted amongst the first passage trajectories to Sf. The 
probability of reaching Sf at time t, before any of the states {S'//}, starting from the state 
Si at time t = 0, is 

roo 

V{Sf\S,)= dtJ{t,Sf\S,), (34) 

Jo 

where J{t, Sf\Si) is the probability flux into the state Sf at time t. Noting that the prob¬ 
ability of jumping directly from the state Si to any other state S is q{S, Si) = ^ ) ; 

and using the backward reasoning, we get, 

P(S/|S.)=«(S/,Si)+ ^ q(S',S,)V{S,\S'). (35) 

The first term is the probability to go directly to Sf from Si and the second is the probability 
to first go to some state S' and then go to S'/ from there (without passing through any of 
the states S//). In other words, the vector P, whose Tth component is P(S'/|Sj), satisfies 
the equation 


M-P = -V, 


(36) 


where V is a vector with components V* = «(S/, S'*). 

Now, the normalized probability distribution of the First Passage Times into the state 
Sf is given by 

F{t,Sf\Si) = J{t,Sf\Si)/V{Sf\S,). (37) 


Using arguments similar to those leading to Eq. (31), the Mean First Passage time can be 
shown to satisfy the following equation: 


5^(>ft)ss.(P(S/|S')r(S,|S')) = -P(S/|Si). (38) 

S' 

For continuous variables, the corresponding Fokker-Planck equation is 

A{si)dsJ^{sf\si) + ]-B{si)ds2V{sf\si) = 0, (39) 

/ ® 

with the boundary conditions V{sf\sf) = 1 and 'P(s/|s/') = 0 for f ^ f. For the MFPT, 

3l(s.)a„(T(s,|Sj)P(ii^|iii) + iB(s,)a,,(T(iiy|Si)P(s/ii.) = -1, (40) 

with the boundary conditions (T(s/|s//)'P(s/|s//)) = 0 for all /', including / = f. The 
applications of these formal expressions are illustrated below. 
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4-. Kramers’method. 


Another method, originally used by Kramers in the famous 1940 paper [SH], can be used 
for the calculation of the Mean First Passage Times and probabilities. In order to calculate 
the MFPT from a state Si to a state 5/, Kramers considered the auxiliary problem with 
an absorbing condition at S'/ and a constant flux J entering at the state S'*. The mean 
time, Tk, that the particles spend in the system, traveling from the state Si to the state S/ 
can be calculated from the average occupancies of all states, N{S,t), which obey the same 
Master Equation as the probability distributions of the individual particles, Eq. Q with the 
extra flux term. Intuitively, in steady state, the flux through the system obeys the following 
relation |49j : 


J=J2 N{S)ITk. (41) 

S^Sf 

It can be rigorously shown that the Kramers’ time Tk is identical to the actual MFPT from 
Si to S'/, proven by by Reimann, Schmid and Hanngi in 153 (see also H). 

Formal derivation. At steady state, the vector of occupancies, N, satishes the equation 
dtN = MfN + J = 0, where Js = JSs,Sr Thus, N{S) = -J{Mj^)sSi and 


Tx = 5^1V(S)/J=-5^(A47‘: 


SSi 


(42) 


from Eq. 41 On the other hand, the vector of the MFPT’s, T, is the solution of the Back¬ 
ward Master Equation Eq.([l 8 |, Ait ■ T = —where I is the unity vector with components 
J 5 = 1 for all S. In other words. 


T(S/Si) = - = Tk, (43) 

s s 

where we have used the fact that the forward operator Ad /■ is the transpose of the backward 
operator: (Mf)ss' = (Mb)s's- 

Kramers’ method can be also extended to calculation of the probabilities, but not the 
times of exit into multiple absorbing states. For instance, the probability to exit through 
state Sf starting from state Si is 


P(SfiS,) = J(Sf)/J, 


(44) 
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where J{Sf) is the steady state flux into the state Sf [58]. Although less general and non- 
generalizable to flnding the probability distributions, Kramers’ method is often a useful and 
convenient way to calculate Mean First Passage Times and probabilities. 


III. APPLICATIONS. 

A. Channel transport. 

1. Background. 

Ubiquitous channels and transporters shuttle various materials into and out of the cell, 
as well as between different cellular compartments. Examples include porins in bacteria, nu¬ 
clear pore complex in eukaryotic cells, transport of polypeptides into the endoplasmic retic¬ 
ulum, ion channels and many others. Their functioning provides inspiration for the creation 
of bio-mimetic nano-transporters for technological applications. During the past decade, 
research of transport through biological and bio-mimetic transporters has seen increased 
application of precise and quantitative biophysical techniques that allow the resolution of 
the durations of the single molecule transport events on the single channel level, in parallel 
with the development of the appropriate mathematical analysis tools. Combination of the 
experimental and theoretical work has resulted in the development of a conceptual frame¬ 
work for the explanation of the transport specificity and efficiency of such nanochannels 

ga EH Eg HU SHlEHHZni . 

Mathematically, transport through a channel can be viewed as a First Passage process 
whose starting point is the entrance of the particle into the channel and its final point is 
the particle exit from the channel. Figure [^ illustrates the different representations of the 
channel transport problem, discussed below. 

2. Discrete channel representation: forward Master Equation method 

Transport of a particle through a channel can be viewed as the hopping between discrete 
sites, as illustrated in Fig. [^ The hopping rates can reflect the energetics, the external 
forces, local geometry or any other factors that affect the particle motion in the channel. 
The model itself is much more general than just a description of the channel transport. With 
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FIG. 2. Channel transport representations, a) Schematic illustration of the channel transport, 
b) Discrete representation of the channel as a sequence of discrete sites, c) Particle movement in 
the channel is viewed as a continuous diffusion in an effective potential U(x). The exit probabilities 
at the channel ends can be represented either via a radiation boundary condition or as absorbing 
boundary conditions located at a short distance a (of the order of the particle size) from the channel 
ends (See text). All these models approximate transport as one-dimensional. Nevertheless, the 
FP methods can be extended to take into account the full three-dimensional nature of the channel 
transport HU. 


the appropriate choice of rates, it has been used to describe molecular motors walking on 
a microtubule, DNA polymerase during transcription, RNA transcript moving through the 
ribosome during translation, or a transcription factor search of the binding site on the DNA 

uni E21 HUTS]. 

The particles start at site i = 1 and hop inside the channel between the adjacent sites 
with the rates until they either translocate through the channel, exiting from site 

N with the rate r_>., or exit the the ’’wrong side” - site 1, with the rate [761180] . The 
probability Pi(t) for the particle to be at site i at time t then obeys the (forward) Master 
Equation 


dtPi{t) = ri_i^iPi_i{t) + ri+i^iPi+i{t) - -f ri^i-i)Pi{t) foYl<i<N 

dtPi{t) = r2^iP2{t) - r^Pi{t) , dtPN{t) = VN^i^^PN-iit) - r^PNit). (45) 
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In the matrix form, 


a,p(t) = M. p(«), 


(46) 


where the vector P(t) = {Pi{t),Pi{t),Pi\f{t)) and the tri-diagonal matrix M has the 
following elements: Mjj = for 1 < j < iV, Mi^i = — (r^ -|- 

^1^2)1 Mn,n = + ^n^n-i)- 

For a particle starting at site i = 1, the initial condition is -Pi(O) = and the solution 

the probability 


to Eq. (46) is Pi{t) = |^e^*P(0)j = . According to Section 

to trans’ 


IID3 


ocate through the channel, exiting through site N, is the integral o 


the probability 


flux out of site N: 


= 




dt = r. { M 


1,N 


r-l 


1,V 


(47) 


The probability density of the transport times distribution is then 




and the Mean First Passage Time is 


1 M~l 

T^ = — tF{t)dt = 

Jo 


(48) 


Special case: uniform and symmetric channel. For a uniform and symmetric channel, 
where all the internal rates are equal, = r, and the exit rates at the ends are equal 

to each other, = r<_ = Tq, the transport probability and the time can be calculated 
analytically |80j : 


- 2 + (N-l)n/r T^^^^(6 + 6Nr„/r + {Nr,/rn 


(49) 


This equation has interesting physical consequences. In the diffusion dominated regime, 
Nro/r S> 1, the probability of translocating is small: ~ nJ/ r because most of 

the particles exit from site 1 soon after the entry, without translocating. In this regime, the 
transport time displays the familiar scaling with the channel length: T_^ ~ N‘^/r. 

By contrast, in the opposite regime, Nro/r -C 1, which corresponds to trapping the 
particle in the channel, the rate-limiting step is the exit from the channel end. In this case, 
the transport time scales linearly with the channel length: ~ N/ro, illustrating the often 
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non-intuitive behavior of the First Passage Times. Despite the fact that the transport time 
is long in this limit, the transport probability increases to = 1/2 independent of the 
parameters. This counter-intuitive fact was hrst realized in the context of the facilitation of 
oxygen transport in tissue by myoglobin [HI]. More recently, facilitation of channel transport 
by molecular trapping, corresponding to small ro/r, has emerged as the explanation of the 
specihcity of channel transport (see also the next section) [IHl |58l [TOl |82l 183] . 

The total mean residence time in the channel, averaged over both translocating and 
returning particles, is Ttot = Note that it scales linearly with the 

channel length, counter to our intuition about the diffusion times. 

In principle, the Mean First Passage Times can be calculated explicitly for any set of 


transition rates either using Eq. (46) and calculating the probability flux, or by solving the 
Backward Master Equation. The hnal answer is obtained in terms of large combinations of 
the transition rates and is very cumbersome. These transport times and probabilities can 
also be obtained using the Kramers method (Section IID4 and [HHl [Hl|)- The methods of 
this section can also be used for the calculation of FPT distributions [HHUHHIHS]- The reader 
is referred to [U ES] EZ] for details; see also Section III B 1 below. 


3. Continuous coordinate representation: backward Fokker-Planck approach. 

Particle motion in the channel can also be represented as continuous diffusion in a poten¬ 
tial U{x) with the diffusion coefficient D{x), which, in principle, can be spatially dependent. 
The discrete and the continuous models can be connected by relating the hopping rates 
between adjacent sites, ri^i±i, to the energy differences: ri±i^i = 

where d is the inter-site distance. However, any choice of rates that satishes the detailed 
balance condition, rj^i/ri^j = is physically acceptable. In the continuous 

representation, the exit rates from the channel at x = 0 and x = L can be taken into ac¬ 
count using the radiative boundary conditions at the channel ends: dx{x,t)\Q = ^p(0,f) 
and dxp{x,t)\L = —-^p{L,t) [511E2]- The constants determine the probability of 

the actually exiting the channel once it reaches the boundary, or getting “reflected” back 
inside. A completely absorbing boundary corresponds to /c = 0, while k = oo corresponds 
to a completely reflective boundary. Thus, they can be related to the rates, r^ and r<_, of 
the discrete case that also reflect the probabilities of the particle at the exit site to leave the 
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channel, 


r+rt- 


and 


r+r_ 


According to Section IID 3, the translocation probability (x) to exit through x = 
L, starting from an arbitrary position x, satishes the stationary Backward Fokker-Planck 


equation (39), compactly written as 

A 

dx 


j = 0, (50) 

with the boundary conditions dxV^{x)\Q = ^'P_s>(0) and dxV^{x)\L = —^V^{L) 


The directional transport times T^{x) can be calculated from the corresponding backward 


equation (40), 


A 

dx 




with the boundary conditions dx{T^{x)V^{x))\o = -j^{T^{0)V^{0)) and 
dxiT^ix)V^{x))\L = -^(T^(L)P^(L)). 

For and 17(0) = U{L), the above equations give for the transport probability, 

F_^iP_( 0 ), 


V^ = 


1 + fc /q /D{y) 

2 + kJ^dyeUiy)/kT/D{yy 

and for the transport time T_> = T_,.(0), 

oL r 


(52) 


k 


1 + 


rx ^U(y)lkT 


D{y) 


dy 


1 + 


L ^U{y)lkT 


D{y) 


dy 


(,-U{x)lkT^^ 


(53) 


Special case: uniform channel. For a uniform potential prohle and constant diffusion 
coefficient, U{x) = E and D{x) = D for all x, one gets for the transport probability 

1 


V^ = 


and time. 


2 + 


L __E6 + 6fe^ + (f)V^ 


Qk^ 


2 + §e^ 


(54) 


(55) 


With the appropriate identification of k, these expressions become identical to the discrete 


channel model, Eq. (49). In the context of channel transport, one is typically interested in 
molecular trapping inside the channels, F < 0. In particular, in the limit of short channel 
and strong trapping, Lke^/D 1, the translocation time ~ Le~^fk is proportional to 
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the channel length and exponentially increases with the trapping energy \E\. Conversely, 
the limit of long channels, Lke^/De S> 1, the transport is dominated by diffusion and the 
transport time obeys the familiar scaling with the channel length L, ~ /D . 

Physical choices of the exit rates and the radiative constants: equivalence of different 
models. The choice of the exit rates in the discrete site method and the fc’s in the radiation 
boundary method depends on the physical problem under consideration. For channel trans¬ 
port, they can be determined from the coupling of the quasi one-dimensional diffusion inside 
the channel to the three dimensional diffusion outside. This can be performed either in the 
forward na or the backward [HS] formalism and results ink = in the radiation boundary 
condition method and rofr = %-^e^ in the discrete site method; a is the channel radius 
and Do is the diffusion coefficient outside the channel. Finally, identifying r = 2D/df {d is 
the inter-site distance), the expressions for the transport probabilities and times obtained 
by the discrete and the continuous methods become equivalent up to a numerical factor of 


4 /7r in the denominator (Eq. (49)). 


4- Mapping onto one-dimensional diffusion. 


Another rendering of the channel transport, which approximates the transport also out¬ 
side the channel as one-dimensional diffusion, is useful for the analysis of transport events 
through individual pores on the single molecule level [HI HTJ [58l [Ml IM!- In this repre¬ 
sentation, the particle starts from the position x = 0 (channel entrance) and performs 
one-dimensional diffusion in the potential U{x) until it reaches an absorbing boundary at 
either x = —a or x = L a, corresponding to the exit from the channel. 


As discussed in Section IID 3, the transport probability (x) to reach L-|-a starting from 


X satishes the stationary Backward Fokker-Planck equation with the boundary conditions 
V^{—a) = 0 and V^{L -I- a) = 1 


OX 




(56) 


which gives 


VPx) = 




( 57 ) 


Note that V^{x) is independent of the diffusion coefficient D. 
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For a flat potential profile, U (x )/kT = E for 0 < x < L (inside the channel) and i7(x) = 0 
outside the channel, assuming that the diffusion coefficient is the same inside and outside 
the channel, the transport probability = '^’->>(0) becomes 




1 

2 + keE/kT- 
a 


(58) 


Note that it is equivalent to the expression obtained in the discrete channel representation. 
The mean transport time obeys the corresponding backward Fokker-Planck equation 

= -V^{x). (59) 

with the boundary conditions V^{L +a)T^{L +a) = V^{—a)T^{—a) = 0. For the negative 
and flat potential prohle, U{x)/kT = E, it yields 


r., = ^ - a/L) + a/L) ~ for |B|AT » 1. (60) 

This expression for the forward time is qualitatively similar to the expressions obtained 
using the discrete and the radiative boundary conditions methods. These results are also 
closely related to the transport of the long chains, such as flexible macromolecules, through 
small pores (see below). 


5. Multiple particles in the channel. 

Until now, we have considered a single particle in the channel. However, a channel can 
contain several particles simultaneously, which interfere with each other’s movement. De¬ 
scription of the movement of an individual particle within the flux of other particles (known 
as the “tracer” particle) is a complicated problem because the motions of the neighboring 
particles are correlated. In this case, there is no closed Master Equation for the probabil¬ 
ity distribution of the tracer particle, and understanding single molecule transport in this 
regime remains a major challenge. 

Exact solutions. It is possible to obtain some exact results for mean residence times even 
for channels with large numbers of particles although the results are typically cumbersome 
[90119^ . Here, we briefly sketch the main points of the derivation for the case of single 
hie transport in a uniform channel in equilibrium with a solution of particles |90]. Most 
generally, the system of multiple particles in a channel is described by the multi-particle 
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probability function P{x,t\y) that the vector of particles’ positions is x at time t, starting 
from the initial vector y [HHl UHl ES] • The crucial insight is that because the particles cannot 
bypass each other, the initial order of the particles is conserved: if ym < yn for sjiy two 
particles at the initial time, it implies that Xm < Xn for all future times. That is, the parts 
of the phase space accessible to these particles are bounded by the planes dehned by the 
condition Xn = Xm in the vector space x. This implies a reflective boundary condition at 
the Xm = Xn plane for any two different particles m and n. 


d„P(x,t\y) = 0 and thus (,d^„P(x,t\y) - d^^P(x,t\y)) = 0 


(61) 


where dn denotes derivative normal to the plane dehned by Xn = Xm- One can then use the 
multi-dimensional generalization of the image method to compute P{x,t\y) and the survival 
probability of the “tracer” particle by integrating out all other coordinates [5l 1901 El] • For 
single-hle transport, the “tracer” particle is known to perform anomalous diffusion with 
the mean square displacement varying with time as (Ax^) ~ instead of the familiar 
diffusion law (Ax^) ~ t. This type of motion can be treated within the anomalous diffusion 
formalism, which, in principle allows calculation of the appropriate First Passage times and 
probabilities [5^ E5] ES] EB]. 

Mean field approximations. Insights into the hrst passage times of interacting particles 
in crowded channels can be obtained using the mean held/effective medium approach that 
approximates the effect of the other particles on the ’’tracer” particle by the average steady 


state density (see [MIIZSIEZIEH]). For the discrete hopping model of Section [III A2| in the 
mean held approximaton the problem reduces to the single particle case with appropriately 
modihed hopping rates, Vi^j —>• — fij), where hj is the average steady state occupancy 

of site i [SUl IM] , Although this method neglects the correlations between the particles and 
misses many important properties of crowded diffusion, it gives reasonable approximations 
for the MFTP. 

For a uniform and symmetric channel with a steady state hux J of particles impinging 
at the channel entrance, the dynamics of the “tracer” particle is then described by the 


discrete random walk model dehned in Eq. (46), with the transition matrix Mi±i^i = r(l — 


hj±i) Mi^i = —r(2 — fij+i — fij-i) and 0 otherwise [80]. The resulting analytical expressions 
are cumbersome, and the outcomes are summarized in Fig. 

The main conclusion is that the crowding increases the average translocation time, while 
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FIG. 3. First Passage Time of the “tracer” particle within steady state flux. Mean 
translocation time of an individual particle within a non-equilibrium steady state flux through 
the channel, normalized by the transport time in an empty channel, T^, as a function of the flux 
through the channel. The lines are analytical results; dots are the simulations. Based on [SOj . 


decreasing the average time of abortive transport events, in which the particle returns from 
site 1. Surprisingly, for the uniform and symmetric process, the transport probability 
and the overall residence time T of the “tracer” particle are the same as in the single-particle 


case of Section III A 2 


= 


T = 




(62) 


2 -h ro{N - l)/r ’ 

This is a consequence of the cancellation of correlations for certain averaged quantities in 
interacting random walks on isotropic lattices Ea. 

Dense regime. The situation simplihes again in the limit of very high densities, when 
essentially all available space is occupied by the particles. This occurs, for instance, for the 
transport of water through nanochannels, such as aquaporins or nanotubes [60l |99l 1100] . In 
this case, the particles can enter and exit the channel only through large collective motions 
of the whole train of particles occupying the channel, whereby the lead particle leaves the 
channel concurrently with the entrance of a new particle from the rear. We denote the 
probability density of the time intervals r between such collective motions as 

The channel can be modeled as a chain of M sites, each occupied by one particle at 
all times. The probability of a particle to be at a position m along the channel at time t, 
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starting from mo, obeys the following Master Equation, 

P(m,f|mo) = 5m,mo(l - / dT'll){T)) 

Jo 

1 


+ - / (ir'0(r)[P(m,f - r|mo + 1) + P(m,t - r|mo - 1)], 


(63) 


with the initial condition P(m, 0|mo) = dm,mo ci-nd the boundary conditions 
P(0,f|mo) = P(M + l,f|mo) = 0, corresponding to the particle exiting the channel. If the 
time intervals between large scale motions obey Poisson statistics with the mean inter-event 
time 1/k, = ke~^^, it can be shown using Laplace Transform that the above equation 

reduces to the familiar random walk on the discrete lattice USD]: 

k 

5tP(m,f|mo) = -[P(m, f|mo + 1) + P(m,t|mo — 1) — 2P(m,f|mo)]. (64) 


This allows the calculation of the transport times and probabilities using the methods de¬ 


scribed in Sections III A 2 and III A 4 [100] : 


1 y M(M + 3) 
M +1 ’ ^ 3fc 


(65) 


Note that the transport probability is again the same as in the non-interacting particle case 
but the translocation time scales as The model also allows to calculate the times of 
more complicated collective motions, such as the interval between the exit times of the hrst 
and the last molecule of the train. This simple model of collective excitations in a strongly 
interacting system is in a very good agreement with the atomistic simulations [100]. 


6. Translocation of long chains through channels. 

The First Passage problem also arises in the context of translocation of long chains, such 
as DNA, RNA and unfolded proteins - and polymers in general - through nanopores. The 
biological examples include translocation of unfolded proteins into the periplasm in bacteria 
and endoplasmic reticulum in eukaryotes. Research on the subject has been driven by the 
technological promise of such devices for DNA and RNA sequencing and protein sorting 

[lai El EH EnMnn. 

If the length L of the polymer is much larger than the thickness of the pore, its motion 
can be viewed as the diffusion of the pore along the polymer, starting from x = L, not unlike 
the models of channel transport illustrated in Fig. Once the pore coordinate x reaches 
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zero, the polymer is considered to have translocated throngh the pore. In the simplest case, 
one can neglect the conhgnrational entropy of the polymer ontside the pore [Ml llU5j . Then 
the probability density of the pore being at a position Q < x < L along the polymer can be 
described by the Forward Fokker-Planck eqnation, 

dtp{x,t) = D{dlp{x,t) - -/-d:,p{x,t), (66) 


with the bonndary conditions p{0,t) = p{L,t) = 0, corresponding to the translocation and 
the retnrn of the chain, respectively; / < 0 is the external force (for instance, electric held) 


that pnlls the polymer throngh the pore m- The general solntion of Eq. ( [6^ for the initial 
condition a; = a;o is 


= I 

n=—oo 


e sm{knx)sm{knXo), 


(67) 


where k^ = ^ and Wn = D [k^^ + 7 (^^) j [51 I1U6] . Using the Poisson identity, Yl'^=-oo /(’^) ^ 
X]m=-oowhere /(27rm) = / dn/(?7,)e*^’^”^"', Eq. (67) becomes 

1 


p{x,t) = - 

which can be rewritten as 

/ N 1 

p[x,t) = 


{x-XQ)f Dt / f \2 
g 2kT g 4 '^kT' 


E 


(re — XQ-\-2Lm) (x4-XQ-\-2Lm) 

g 4Dt — g 4Dt 


( 68 ) 


yjA'xDt _ 


E 


fL 

ekT' 


— {x — XQ—^^t-\-2Lm)^ 


(x + xo-§ft+2Lmp fm 

— e 


(69) 


Each term in these inhnite series can be interpreted as an “image” particle with the starting 
point at —Xq, 2L — xq, —2L + Xq, 2L + xq etc., snmmed with the appropriate weights as 
to satisfy the bonndary conditions p(0, t) = p{L^ t) = 0, analogons to the solntions to the 
Poisson eqnation in electrostatics jSlEaiM]. 

The translocation probability can be calcnlated exactly: 

,.00 . _ f(L-^o) 


= 


J(0, t)dt = 


fL 

1 — e kT 


(70) 


where the probability hnx into the absorbing bonndary at rc = 0 is J{0,t) = \Ddxp{x,t)\x=o 


(see Section IID). The normalized probability distribntion of the translocation times, F{t), 
and the average translocation time, are 


F{t) = J{0,t)/V^ , = lim 

xq^L 


tF(t) , 


(71) 
























which result in rather cumbersome expressions. However, the probability distribution of the 
translocation times can be approximated (for LF‘/[Dt) S> 1) as 

2 


F{t) ~ 




{DtfF \Dt 


/{ADt) 


(72) 


which has a maximum around fmax = ^^(1 ~ ^17^ + ...) |in6j . The maximum of the 
probability density is an alternative characteristic of the typical translocation time. Note 
that for heavily asymmetric distributions, it can differ significantly from the mean time. 

Approximations. For strong forces, or long channels, the typical translocation can be 
viewed as an almost deterministic motion in the direction of the force, with the mean 
“velocity” v = In this case, the probability that the chain does not translocate is low, 
and one can move the absorbing boundary condition at x = L to cx) |l2l [SH 110811109] . This 
greatly simplifies the problem, which now has only one absorbing boundary condition at 
X = 0. Taking the limit L —)■ cx), Eqs. ( |^ and ( [M| ) reduce to [SllM] 


p{x,t) = 


3 kT j 


y/AirDt 

The distribution of the First Passage times is 

F{t) = lim J(0, t) = 


(^+^0-^) ^0/ 
g ADt g kT 


(73) 




xq — vL 


ty/AirDt 


(74) 


and mean translocation time is 


= 


TF{T)dT = 


kT l + 


(75) 


D \f\Ll-eFf\L/kT 

As expected, for strong bias, 
|/|L/fcT 3> 1, the translocation becomes an essentially deterministic motion with velocity 
V = D\f\/kT, so that T_> ~ L/v ~ fmax [3 11011 IllOj . 


[3- Note that to the first order in it is identical to 


B. Receptor binding and adhesion. 

Another class of phenomena that are naturally described in the First Passage process 
formulation is the multivalent binding and adhesion - from macromolecular association to 
receptor signaling and viral cell entry [HI 111111116] . In this section we review several recent 
works illustrating the applications of the First Passage methods to these problems. 
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1. Viral particle binding and dissociation at the cell surface. 


Typically, the first stage of viral entry into a target cell is the binding to the cell surface 
receptors. The lifetime of a virus particle (virion) on the surface of a target cell is an 
important early determinant of the infection outcome. In the model of |117j . the virion has 
N sites on its surface that can bind receptors on the cell surface; the latter are present in 
the surface concentration C. The virion is thus in one of the N states: with n = 1, 2, 3, ...N 
sites bound. The state with n out of N sites bound can transition into the state with 
n — 1 bound sites, through breaking of one bond, with the rate fco for n = 1 and nk^i for 
n > 1. Alternatively, any of the unbound sites can form a new bond with a surface receptor, 
resulting in a transition to the n + 1 state, with the rate {N — n)kiC. Physically, the 
rates ki and k_i reflect the local “on” and “off” rates, primarily determined by the binding 
energy, while fco reflects not only the the time of the local bond breaking but also the time of 
diffusing away from the cell surface. The process is illustrated in Fig. Note the similarity 
with the kinetic scheme for the particle in the channel of Fig. Similar models have been 
used to describe nanoparticle adhesion onto cell surface |118] . 


In the backward approach of Section IIC, the mean time to unbinding starting from n 
bound sites, T„, satishes the equation 

1 A, 


T = 

n. 


-Tji-i + 


/^7i 


- n+1? 


(76) 


T A^ A^ T /ijj A^ T /ij, 

where Xn = {N — n)kiC, fin = nfc_i, fii = ko |11711119] . 

Typically, the virion binding starts from just one site, n = 1. In this case, the sequence 


of difference equations (|76|) can be solved analytically, giving for the average lifetime of the 

■(l + i^C')^-l' 


virion on the cell surface 


Ti = h 

ko 


NKC 


(77) 


where K = ki/k^i is the affinity of an individual site to a surface receptor 1119] . 

Physically, the binding affinity K is the inverse of the dissociation constant and is 
related to the binding energy, e > 0, as iP ~ the dissociation constant Kd is sometimes 

colloquially referred to as the “affinity” as well. In the limit of weak binding or low surface 
receptor density, KC <C 1, the dissociation time is Ti ~ l/ko, indicating that the virion 
is most likely to escape immediately after binding without recruiting additional surface 
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receptors. In the opposite limit of the very strong binding, S> 1, Ti ~ 

and is exponential in the binding energy and the number of the binding sites, indicating 

that in this limit it is essentially cooperative binding that engages all the binding sites 

simultaneously. 


2. Multivalent binding: avidity 


Similar problems arise in many instances of multivalent binding. For example, trans¬ 
port proteins shuttling cargoes through the Nuclear Pore Complex (NPC) possess multiple 
binding sites to the hydrophobic residues located on the natively unfolded proteins located 
within the Nuclear Pore Complex [emiMi]. First Passage theory can be used to analyze the 
results of single molecule fluorescence tracking experiments to infer the binding times and 
the the effective affinity of the transport factors to the hydrophobic repeats |11]. Assuming 
that the hydrophobic repeats are present in volume concentration F within the lumen of 
the NPC, the problem of the transport factor binding to the NPC becomes mathematically 


identical to the previous section. The dissociation time can be calculated from Eq. (77), 


which defines the effective “off” rate of the interaction. Together with an “on” rate of the 
first binding event, kon, the effective dissociation rate determines the effective interaction 
affinity: iPeS = NkonTo^i- The factor N arises because any one of the N binding sites on the 
transport factors can bind a hydrophobic repeat first. 


For a transport protein with four binding sites, expanding Eq. (77), 


A'.j = Ko(iF-KF+ (KFf + 


(78) 


where Kq = kon/ko is the affinity of the hrst binding event [H]. This effective affinity 
Aeff cannot be derived solely from the binding energies, but depends on the concentration 
and the availability of the binding factors - effect known as the “avidity” for multi-valent 
interactions [nn na nni EH. 


3. Competition between viral dissociation, endocytosis and fusion. 

Other possible outcomes of the virion binding to the cell surface, in addition to dissocia¬ 
tion, were considered in |113] . While the virion is bound to the surface receptors, it can fuse 
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with the cell membrane and deliver its genetic material into the cell. On the other hand, it 
can become engulfed by the cell membrane, endocytosed and targeted for destruction. The 
fates of the virus and of the infected cell are determined by which of these three processes 
completes hrst. The kinetic scheme of the model is illustrated in Fig. The virion can be 


^ Unbidning Fusion {N-1)k^ Fusion k Endocytos 

O 'Qi (Q) 


FIG. 4. Virus binding to the cell surface. Schematic representation of the the kinetics of virus 


binding to the cell surface, reviewed in Sections III B 3 and III B 1 In the model of Section III B 1 

ke = kf = 0. 


in one of the N states with n bound receptors plus the non-specihcally adsorbed state at 
n = 0, from which it can completely dissociate with rate kd- In addition to the transitions 
from the state with n bound sites to a state with u ± 1 bound sites and the dissociation from 
state n = 0, the process can terminate from n = N through endocytosis with rate or 
via membrane fusion from any state n with the rate nkf. One important difference of this 


model from Section |III B 1[ is that only binding sites that lie close to the circumference of 
the bound area of the virion can bind or unbind. Thus, the transition rates from state n to 
n — 1 (unbinding of one receptor) and to n + 1 (binding an additional receptor), Qn and pn, 
respectively, are Pn/ki = ~ ) = (|^(1 - (1 - 2 n/iV) 2 )V 2 for iV > 1. 

The Pn and qn correspond to the and the from Section III B 1 Overall, there is a race 
between the three possible outcomes: dissociation, endocytosis and fusion that occur with 
the corresponding probabilities, that sum to one: = 1. 

The Forward Master Equation for the probability to be in state n, Pn{t) is 


dtPnit) = -{nkf +Pn + qn)Pn{t) + qn+lPn+l{t) + PnPn-l{t) for 0 <n < N 

dtPo{t) =—{kd + Po) + QiPiit) ( 79 ) 

dtPN{t) = -{ke + qN)PN{t) +PN-lPN-l{t)- 
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Using the numerical solutions of the Forward Master Equation, the probabilities of different 
outcomes can be calculated using V-^ = kfU Pn{t)dt, PN{t)dt and = 


kaj^Poit) ms]. 


These probabilities can also be calculated directly using the Backward Master Equation 


method. For instance, following Section IID3, the fusion probability obeys the following 
set of equations: 


-nkf = -{nkf +Pn + qn)Pl + PnVl+i + qnVi-i for u > 1 

0 = -(fed + po)Vl + pqVI (80) 

-Nkf = -{Nkf + ke + qN)pN + 


which can be solved by any method for solution of systems of linear algebraic equations. 

Even this simplihed model predicts rich behavior with important biological implications. 
It be extended to include co-receptor binding and viral exocytosys from the infected cells 

[851 [T^. 


C. Single-cell growth and division. 

1. Biological context. 

The interplay between lengthscales and timescales in the context of cell growth and 
division can be cast as a First Passage Time problem by formulating how interdivision 
times are informed by the stochastic increase in cell size of individual cells. While the 
questions and modeling challenges in this context have been long appreciated |123H135] . 
there is renewed interest [13611164] due to the recent availability of large datasets for single 
cell growth trajectories and cell divisions, made possible by major breakthroughs in single 
cell technologies for unicellular organisms [HOllTT] . 


2. Formulating cell division as a first passage time problem. 

During each interdivision period, the size of a cell (assumed proportional to its mass 
[1251ll65[ll66j j increases according to “the growth law” on average |124[ 112511167] . Typically, 
this increase is either linear or exponential for unicellular organisms [30l [3U 11251116711169] . 
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There are three commonly considered scenarios for how the stochastically increasing cell 
size could inform the cell division. They are [301 11251 llf>71 11681 I17U] . (i) “absolute size 
thresholding” (the “sizer model”), in which the (stochastic) cell size attains a critical or 
threshold value at division; (ii) “differential size thresholding” (the “adder model”) in which 
the thresholded variable is the change in cell size from its initial to hnal value; and, (hi) 
“ratio size thresholding” (the “timer model”) in which the ratio of the size at division to the 
initial cell size is thresholded. For clarity in elucidating the methodology, here we assume 
that the cells are in balanced (steady state) growth conditions, and that intergenerational 
correlations are negligible. Together they imply that the statistics of growth and division 
are identical and independent for all generations of the cells. 

Under these assumptions, the formulation of cell division as a hrst passage time problem 
requires that the following be specihed: (i) a stochastic model for how cell size, s, increases 
with time, t, between divisions; (ii) the function of the cell size, s, that attains critical 
or threshold value, 6, at division; and (iii) appropriate initial conditions, including the 
initial distribution of cell sizes. While in all previous cases considered, every member of the 
ensemble (i.e., each cell in the population) was assumed to be subjected to identical initial 
conditions, in this section we relax that condition to allow different cells to experience 
different initial conditions (such as different freshly divided cells having different sizes). 
This is an added source of stochasticity, namely, extrinsic noise in addition to the intrinsic 
fluctuations encoded in stochastic growth for a given initial condition. 

The time, t, is thus equal to 0 for a newly divided cell and to r at hrst passage, i.e., at 
division. The goal is to then compute the hrst passage time distribution, F{t), where r is 
the interdivision time, i.e., the time taken for the thresholded variable (cell size or function 
thereof) to reach the threshold value, 6. In this section, for additional clarity, we explicitly 
write out the parameters on which the FPT depends in its argument, separated from the 
variable, r, by semi-colons. Thus, division at a threshold 6 has the FPT distribution: F{t] 6). 


3. Relation between cell sizes at division and interdivision times. 

A convenient simplihcation of the FPT problem is obtained using a generic feature of 
cell growth: cell sizes always increase monotonically with time (for living cells), even though 
the increase is stochastic [30l I138j . The possibility of multiple crossings of the threshold is 
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then automatically ruled out (since the threshold is passaged exactly once) and so we do 
not need to consider the auxiliary problem with absorbing boundary conditions (see Section 


IID). This leads to the following mathematical identity. Quite generally, when a stochastic 


variable s increases monotonically with time, its time dependent distribution, P{s, t), can be 
related to the distribution of First Passage Times, F{t] s = 6), through a simple geometric 
argument, which illustrates the derivations of Section |II D[ From Fig. using probability 
conservation, it follows that the cumulative of the size distribution at the threshold value 
must be equal to the cumulative of the First Passage Time distribution. Thus the FPT 
distribution can computed using the relation: 


F(r; 9) = dr 


ds P{s, t) 


(81) 


Je 

When a discrete growth model is used for s, the integral should be replaced by an appropriate 
sum. In this section we denote cell size by s, irrespective of whether the stochastic growth 
model used is discrete or continuous. 



FIG. 5. Cell division as a First Passage Time problem, (a) Schematic of stochastic cell size 
increase from a common initial condition. Between times r and r+Ar, the blue growth tracks cross 
the threshold size, 9. The green trajectories cross the threshold before this interval, and the red 
trajectories cross the threshold after this interval. Using probability conservation, the cumulative 
probability that the size is greater than 9 (above the black dotted horizontal line) must be equal to 
the complement of the cumulative probability that the First Passage Time is less than or equal to 
T (left of blue dotted vertical line at r). (b) Scaling of the First Passage Time distribution. The 
shape of the mean-rescaled division time distribution is timescale invariant, i.e., independent of k, 
when there is a single timescale, 1 /k oc (r), in the FPT dynamics. 
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4-- Scale invariance of the FPT distribution. 


Division time distributions from different growth conditions, for the same organism, have 
been observed to undergo scaling collapses, when rescaled by their condition-specific mean 
values [301 11391 1171] . This observation encodes a deeper truth about FPT distributions: 
whenever a single timescale dominates the stochastic growth and division dynamics, irrespec¬ 
tive of the functional from of the growth law, or the thresholding scheme, the mean-rescaled 
FPT distribution from different growth conditions is scale invariant. 

To see how this result arises formally, we denote the assumed single timescale in the 
problem by k~^ (in practice, this timescale can be tuned by external parameters). For the 
growth variable, s, whose time dependent distribution is P{s,t), and 6 the threshold at 
which division (first passage) occurs, the FPT distribution is given by: 


F{t; 9) = dr 


^ d(^n t) 


ds P{s, r; n) 

poo 

/ ds P(s, kt; K = 1) 


(82) 


Thus, if we now change variables to f = kt and look at its probability distribution, F(r; 9), 
then we have, 

F(fi«) = t 

K 


F{T-9) = df 


ds P{s, r; K = 1) 


Ue 


(83) 


which is manifestly k independent. Therefore, F{f]9) is timescale-invariant and is the 
functional form of the scaling invariant mean-rescaled FPT distribution. 

It is also straightforward to show that the mean FPT, (r), is proportional to 1/k. Thus 
the mean-rescaled division time distribution from different growth conditions will be found 
to undergo a scaling collapse, provided the underlying stochastic growth model has only one 
timescale (as is true for exponential or linear growth), and the thresholding does not itself 
introduce new timescales into the FPT process (as assumed above). 

Conversely, the observation of a scaling collapse of mean-rescaled FPT distributions con¬ 
firms that a single timescale governs the underlying stochastic dynamics. While the formal 
result appears to be intuitively obvious, the implications of observing this in a real biological 
system are significant: the growth law, regardless of its functional form, must depend on 
just one timescale; the thresholding scheme does not itself introduce a new timescale into 
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the division dynamics; the division and growth timescales must therefore be proportional to 
each other and the mean division time, as external parameters are changed. 

The specihc stochastic growth models that we consider in the following sections all have 
one timescale governing growth, and the thresholding schemes we have enumerated above 
are timescale independent. Thus, the division time distributions for each case are scale- 
invariant, when mean-rescaled. (This can be checked directly from the analytical forms 
derived below.) 


5. Master Equation approach. 


In this section we hnd the FPT distribution using the Master Equation approach, by 
using discrete stochastic growth models for s and the Master Equation framework (see 


Section IIBl). In the continuum limit of s, the results for the FPT distributions are 


essentially unchanged. As previously noted, for the purpose of the present discussion we 
assume that cell size growth is either linear of exponential. (Also see “phase oscillator 
model” below for a different interpretation of the linear growth model.) 

Linear growth: In a simple stochastic (discrete) model for linear growth, represented by 


^ , -I 

s -A s -|-1, 


(84) 


the time evolution of s is governed by the Master Equation (see Eq. ([T 


dtPis,t) = k[P{s - l,t) - P{s,t)]. 


(85) 


(See, for instance, [2].) We have made the standard assumption of exponentially distributed 
waiting times. Using standard techniques [2], it is straightforward to show then that the 
ensemble mean, {s(t)) = grows linearly with time as kt. 

We now consider different scenarios for thresholding the size variable in this model. The 
initial condition is that the cell size distribution at t = 0 is R{sq), where s = Sq at f = 0. 
We hrst consider the simplest initial condition: all cells start out with the same s = 0 
at f = 0, i.e., R{so) = 5so,o- If can be shown that for this initial condition the time- 
dependent distribution of sizes at any given time is then the Poisson distribution |21 1172j 
whose single (time dependent) parameter equal to both the mean and the variance of the 
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Poisson distribution, is kt: 


P{s,t) 


e {kty 


( 86 ) 


Now, for the absolute size threshold, s = 6, the FPT distribution is found by using equation 


F(r; 9) = dr 


.s=e 


(87) 


Upon evaluation, we hnd that the FPT obey a Gamma distribution whose shape parameter 
is given by the magnitude of the threshold, 6: 


F{T-9-k) 


ke {kr) 


( 88 ) 


where P[a:] is the Gamma function |173] . and we have explicitly written out the parametric 
dependences of F, as previously mentioned. Note that the FPT distribution F(r; 6] k) = 
kP{9 — 1, r), which is the probability flux from state 0 — 1 to state 9, in accord with Section 

umi 

When all cells are assumed to start with the same initial size, Sq > 0, the size and division 
time distribution for absolute thresholding are, respectively, a shifted Poisson and a shifted 
Gamma distribution, as one might intuit: 


F(s,t;so) 


(s - So)! 


0(s - So), 


(89) 


where 0(s — sq) is the Heaviside Theta function. The corresponding FPT distribution is 


again found using (81): 


F(r; 9 - sq; k) = -—--^-, for 6^ > So + 1 and 

P[6' - soj 

= k for 0 = So + 1, 


(90) 


which is again identical to the probability flux into state 9. A limiting case of the above 
solution, 6* = So + 1, is consistent with the assumption that the waiting time distribution is 
exponential. In this model the initial value, so, does not affect the propensity for stochastic 
growth. Thus, when there is an initial distribution of sizes, denoted by F(so), the result¬ 
ing FPT for absolute size thresholding is given by the convolution of the above Gamma 
distribution with F(so). 
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Next consider differential size thresholding, i.e., a cell with an initial size sq divides when 
its size reaches s = Sq + where the additive threshold, A6 is a given positive number. 
Using the result above for a given initial condition, sq, we hnd now that the division time 
distribution, for additive thresholding: 




ke {kr) 


X 


-1+A0 


ke (kr) 


-i+Ae 


S0=0 

oo 


T[Ae] 


r[A0] 


-, for A9 > 1, 


Riso) X ke-^^ = ke-^A for AO = 1. 

so=0 


(91) 


Not surprisingly, R{so) drops out of the expression, and the FPT distribution for this case 
is thus independent of the initial size distribution. 

In all the above expressions for FPTs, we could set k = 1, i.e., effectively measure all the 
times (including the division time, r), in units of 1/k. Evidently, when this substitution is 
made, the FPT distribution becomes timescale invariant. Moreover, in all cases, (r) (x 1/k. 
Thus a scale-invariant result is obtained when the FPT distributions are rescaled by (r), in 


agreement with the general scaling result derived previously in Section III C 4 


Exponential growth: A simple model in which the ensemble average of s grows exponen¬ 
tially with time, represented by the growth process. 


ks , 

s — y s -|-1, 


undergoes time evolution governed by the Master Equation (see Eq. 0 ). 

dtP{s, t) = k [(s — 1)^(5 — l,t) — sP{s, t)]. 


(92) 


(93) 


As in the previous section, we hrst consider the initial condition where all cells have the same 
initial size, SQ) and then generalize to the case where they may have a distribution, P{so). 
We hnd that the cell size distribution is a Negative-binomial distribution |138j (which is the 
discrete analogue of a Gamma distribution) |174] : 

F(s, r; So; k) = ^ “ M (1 - e-^""°0(s - Sq). (94) 


The FPT for absolute size threshold 9 is then found using (81) to be the Beta-exponential 
distribution |175] . 


F{T;so,9;k) = 


_ g-fcT^(-l+6»-so) 

/?[so, (0 - So)] 


(95) 
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where I3[x,y\ is the Beta-function |173j . As before, it can be shown using the properties 
of the Beta-function that F(r; Sq, 6] k) = kP{6 — 1, r; Sq! k), the probability flux into state 


6*, in accord with Section IID For an alternative approach leading to the Beta-exponential 
solution, see ing. Note that the FPT distribution, which is a Beta-exponential in r, is 
actually a Beta distribution if one transforms the variable t to v = e~^'^. The parameters 
of this Beta distribution are restricted such that the FPT distribution is always unimodal 
in this problem. 

We now formulate the ratio thresholding problem for exponential growth. Starting with 
an initial size, sq, drawn from an initial size distribution, R{so), each cell is assumed to divide 


when its size reaches a ratio threshold, s(r)/so = r. Upon solving the Master Equation (93) 
we get: 


P{s,t;k) = J2 Ri^o) 

so=0 


X 


s - 
■So 


(1 _ 


(96) 


Therefore, using Eq. (81) again, the FPT using a ratio threshold, r, is 

_ g-A:r^(-l-|-rso-so) 


F{T;r-,k) = ^ R{so) 

so=0 


X 


^[so, (rso - So)] 


(97) 


whose shape depends on the specihc choice of initial size distribution, R{so). 

We note that it may be possible to invoke overarching biophysical principles constraining 
growth and division in population balance to self-consistently determine the initial size 
distribution, R{sq). However, this discussion is outside the scope of this review. 

To recapitulate, we have shown in this subsection how given a (mean) growth law and a 
thresholding scheme for division, the FPT problem for cell division can be formulated and 
the cell division time distribution can be analytically derived. Such a model may be used 
to make other predictions about the biological system, including placing constraints on the 
possible topologies of the networks governing the stochastic growth |138] . All division time 
distributions derived above, for different growth laws and thresholding schemes, are posi¬ 
tively skewed (have a long right tail), and are unimodal. However, it is worth mentioning that 
even with the kind of high quality data available in recent single-cell experiments laniEi], 
using the shape of the observed division time distribution to infer the underlying growth law 
(for example, to distinguish between linear and exponential growth) is extremely challenging 
and not practically feasible. 
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6. Cell cycle as a phase oscillator. 


In some scenarios the cell cycle (i.e., the intervening period between successive divisions) 
is modeled as a phase oscillator, with the cell cycle phase, 0 being set equal to 0 for a newly 
divided cell and 0 = 27r for a cell about to divide |125| 11771 - 1179] . If N sequential steps need 
to be completed, as the cell cycle phase increases from 0 to 27i in steps of 271/N, and if 
the waiting time distribution for each step is exponentially distributed as kexp{—kt), then 
the FPT problem for cell division in this model is essentially identical to the one for the 
stochastic discrete linear growth model solved previously. Thus the FPT distribution is the 
gamma distribution: 

U p-kr /L 'jTV-l 

F(r; N] k) = --, for iV > 1 and 

= k hr N = 1. (98) 

In principle, the observed cell division time distributions can therefore be htted to a gamma 
distribution, and used to estimate the number of “elementary” steps in the cell cycle, N. 
However, as we have shown above, the model makes the simplistic assumption, almost 
certainly violated by all realistic systems, that all steps have the identical and exponentially 
distributed waiting-time statistics. Thus, the practical utility of such an estimate of “N” is 
limited. 


7. Phenomenological approach using Langevin or Fokker-Planck frameworks. 


In the examples that we have considered thus far, we have used a microscopic model to 
motivate the growth law, and then used it to hnd the FPT distribution. A complementary 
approach is the phenomenological one in which one proposes Langevin dynamics consistent 
with the observed mean growth law, and assumes an ansatz for the noise term for cell size 
fluctuations. By then going to the corresponding Fokker Planck description, one can use 


standard techniques |5j (see Section IID) to compute the FPT distribution. 

We elucidate this alternative approach with the specihc case of the exponential growth 
law, with ratio size thresholding. This methodology can be readily adapted to other growth 
laws and thresholding schemes. 

The cell size increases, on average, as {s{t)) = (s(0)) exp(H); the cell divides at a time. 
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r, such that s(r)/s(0) is a constant. Motivated by the exponential growth law and ratio 
thresholding scheme assumed, we dehne a new stochastic dynamical variable for each cell, 
x(t), as: 


x{t) = log [s(t)/s(0)], 


(99) 


where s(0) is the initial size of the cell under consideration. The threshold for division is 
then given by 


Xo = x{t) = log 


L.(0) 


( 100 ) 


We then write a stochastic growth model as a Langevin equation (see Eq. (14)) for the 
“Brownian motion” of x{t): 

dx(t) 


dt 


= ft + Vb x(t), 


( 101 ) 


where k is the mean (ensemble averaged) exponential growth rate, the “drift” term; the 
second term is the ‘noise’ term: x{t) is standard Gaussian white-noise, and the “diffusion” 
term measures the strength of the noise in x. To evaluate the FPT distribution, we hrst 


recast this Langevin equation to its equivalent Fokker Planck equation, Eq. 12 and Eq. (13), 
using standard techniques. It is: 




( 102 ) 


For computational ease, without loss of generality, we shall assume that each cell (the 
“random walker”) starts out at x = Xq at t = 0 and divides (gets absorbed) when it hrst 
crosses the origin, x = 0. Thus our initial condition is that P{x,t = 0) = 6{x — xq) and 
we shall impose absorbing boundary condition at x = 0. The fraction of random walkers 
disappearing at x = 0 between times r and r -|- dr is then related to the current of walkers 
entering x = 0 in that time interval according to the probability conservation equation. 


Eq. (21), 


F(r) dr = J [P{x,t) — P{x,t + dr)] dx 
= r P{x,T)dx\ dr, 


(103) 
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where F{t) is the First Passage Time distribution sought. We note that P{x,t) is “normal¬ 
ized” at each time such that 


dxP{x,t) = l— / cItF^t), 


(104) 


while the FPT, F{t), is itself correctly normalized to 1. This is because P{x,t) quantihes 


the density of the surviving random walkers (see Sections IID and IID4). 


The solution to the Fokker Planck equation (102), with the specihed boundary condi¬ 
tions, is most elegantly computed using the method of images, which is routinely used for 
electrostatics problems with symmetry [HllSl]. We note that, contrary to naive expectation, 
the random walker and its image must move in the same direction to satisfy the Fokker 
Planck equation. The solution, which can be verihed straightforwardly by substitution, is 

1 


P{x,t) = 


y/27rB t 


^—{x—XQ-\-Kty/2Bt 


(1 


^—2xxolBt\ 


(105) 


Using this expression for P{x,t) and Eq.(103), we hnd the First Passage Time distribution. 


F(r) = 


Xo 


V2ttB ' 


^-{xo-kt)'^I2Bt 


(106) 


which is the so called Inverse Gaussian Distribution, with parameters k and B. It was hrst 
derived as the FPT distribution for Brownian motion by Schrodinger |18Uj . 

Notably, the shape of the Inverse Gaussian Distribution “scales” with the mean, as ex¬ 
pected (see preceding section on scaling of FPT distributions). As previously noted, cell sizes 
are observed to increase monotonically, even when fluctuations are considered. This implies 
that for the phenomenological descriptions such as presented above (which assume Gaussian 
white noise), the drift term must be overwhelmingly larger than the diffusive term, resulting 
in (approximately) monotonic growth. In this drift dominated regime, the Peclet number 
for this FPT problem j5] is thus a very large dimensionless number.For this model, it is equal 
to l/? 7 ^, where r] is the coefficient of variation of F{t). In this regime kcho 3> \/B ao and so 
the “distance” travelled by the peak of the distribution of x in a given time interval is much 
greater than the corresponding widening of the distribution. In the drift dominated limit 
the problem is approximately equivalent to the motion of a Gaussian probability packet, 
Q{x,t), which is the solution to the same problem with the same initial condition, but with 
no absorbing boundary condition at x = 0. The error in this approximate solution comes 
from the fact that we do not take into account that the same particle might have crossed 
the boundary at x = 0 more than once. 
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We note that parameters of this phenomenological description can be inferred from ex¬ 
perimental observations. The drift, k, is directly given by the ensemble-averaged growth 
curves. From time series growth data, the mean squared displacement of x with time can 
be computed to conhrm that the behavior is diffusive, and to read off the diffusion strength, 
B. The scale invariance of the FPT, when external parameters are tuned, provides an ad¬ 
ditional check on whether these two parameters are then found to be related to each other 
as predicted. 


IV. CONCLUDING REMARKS. 


Recent years have seen a renewed interest in stochastic processes in biology, including 
First Passage processes, resulting in a rapidly increasing wealth of literature. The aims of 
this review were two-fold. First, we consolidated the theoretical foundations and techniques 
of FP processes within a unihed framework. Second, we provided an introduction to the 
practical use of FP methods in biophysical applications using several pertinent examples. 

Out of necessity, the applications discussed here do not constitute an exhaustive list 
of potential uses of the FP theory, even just in the cellular context. We apologize to the 
authors whose work could not be cited. However, the methods and techniques discussed here 
have been successfully applied to molecular motors, translation, transcription, protein and 
enzyme dynamics as well as signaling. Beyond the context of cell biology, vast literature on 
FP applications to neurobiological systems and population genetics can be found. We point 
the interested reader to some relevant literature on these topics [IHl ESI ESI [TH ll81H189j . 

Several important theoretical aspects also were not reviewed here, most notably those 
pertaining to systems with fluctuating barriers and boundaries [HHl I19UI1191] , anomalous 
diffusion [521 ESI 1192] and Hamiltonian methods for large fluctuations [521 1186] . Finally, in 
this review we focused entirely on the analytical approaches. Simulation techniques that 
complement analytical approaches have played a crucial role in revealing the stochastic 
mechanisms of cells and molecules. We direct the reader to the following textbooks and 
reviews as starting points for further inquiry [5TI I182[ 1193] . 
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